function fig7_18

global a b

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0
%  where y = (y1, y2, y3 , ..., yn) is an n-vector

clf
% get(gcf)
set(gcf,'Position', [17 938 667 399])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

for i=1:2

    if i==1  % comes close to (0,0)
        tmax=400;
        a=-1; b=1;
        delta=0.1;
        y10=1; y20=0;
        y0=delta*[y10 y20];
    else   % comes close to (0,0)
        a=1; b=1;
        y10=1; y20=0;
        c=1/(y10^2+y20^2);
        delta=0.1;
        tmax=0.99*4*c/(a*delta^2)
        %tmax=390
        y0=delta*[y10 y20];
    end

    %  calculate solution using a MATLAB routine
    opts = odeset('RelTol',1e-12,'AbsTol',1e-12);
    [t,y] = ode45(@rhs,[0 tmax+5],y0,opts);
    %[t,y] = ode23s(@rhs,[0 tmax],y0);
    if i==1
        subplot(2,1,1)
        hold on
        plot(t,y(:,1),'b','LineWidth',1.0)
        c=1/(y10^2+y20^2);
        tt=linspace(0,tmax,200);
        for ii=1:200
            A(ii)=delta/sqrt(c-0.25*a*delta^2*tt(ii));
        end
        plot(tt,A,'r','LineWidth',1.1)
        axis([0 tmax -0.1 0.1])
        grid on
        box on
        set(gca,'ytick',[ -0.1 0 0.1])
        set(gca,'YTickLabel',{'$-0.1$';'0';'0.1'})

        ax = gca;
        ax.XAxis.TickLabelInterpreter = 'latex';
        ax.YAxis.TickLabelInterpreter = 'latex';
        set(gca,'FontSize',13)

        xlabel('$t$\,-\,axis','Interpreter','latex','FontSize',14)
        ylabel('$x(t)$','Interpreter','latex','FontSize',14)
    else
        subplot(2,1,2)
        hold on
        plot(t,y(:,1),'b','LineWidth',1.0)
        c=1/(y10^2+y20^2);
        tt=linspace(0,tmax,200);
        A=zeros(200,1);
        for ii=1:200
            A(ii)=delta/sqrt(c-0.25*a*delta^2*tt(ii));
        end
        plot(tt,A,'r','LineWidth',1.1)
        axis([0 tmax+5 -1 1])
        grid on
        box on
        set(gca,'ytick',[ -1 0 1])
        set(gca,'YTickLabel',{'$-1$';'0';'1'})

        ax = gca;
        ax.XAxis.TickLabelInterpreter = 'latex';
        ax.YAxis.TickLabelInterpreter = 'latex';
        set(gca,'FontSize',13)

        xlabel('$t$\,-\,axis','Interpreter','latex','FontSize',14)
        ylabel('$x(t)$','Interpreter','latex','FontSize',14)
    end
    %plot(y(1,1),y(1,2),'.b','Markersize',20)
end

%exportgraphics(gcf,'/Users/mark/Desktop/stable2.pdf')


%  define f(t,y)
function dy=rhs(t,y)
global a b
dy=zeros(2,1);
dy(1) = y(2)+a*y(1)*y(2)^2;
dy(2) =-y(1)+b*y(1)*y(2)^2;













